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Abstract. In this paper we describe a method to compute Generalized Polarization Tensors. 
These tensors are the coefficients appearing in the multipolar expansion of the steady state 
voltage perturbation caused by an inhomogeneity of constant conductivity. As an alternative 
to the integral equation approach, we propose an approximate semi-algebraic method which is 
easy to implement. This method has been integrated in a Myriapole, a matlab routine with a 
graphical interface which makes such computations available to non-numerical analysts. 



The classical electrical impedance tomography problem in two dimensions can be described 
as follows. Suppose that C K 2 is a bounded simply connected domain. Given a conductivity 
map 7 6 L°°(Q), suppose that there exists a constant c > such that c _1 < 7 < c. Define the 
Dirichlet-to-Neumann map A : J ff 1 / 2 (9fi) -> H- 1 ^ 2 (dD,) by 



where u <E H l (Q) is the voltage potential associated with 7, that is, the unique solution of 



The electrical tomography problem is then to reconstruct 7 from A. The fact that 7 is uniquely 
determined by A has recently been established in dimension two |10| . and this question is still open 
in dimension three. It is known however that, without additional assumptions on the conductivity, 
this problem is extremely ill-posed This ill-posedness can be considerably reduced in several 
cases of practical importance by the introduction of a priori information on the form of the 
conductivity to be recovered. One such problem is the following: given complete (or incomplete) 
knowledge of the Dirichlet-to-Neumann map, estimate the internal conductivity profile under the 
assumption that 7 is constant except at a finite number of small inhomogeneities. This problem, 
and some variants, has been investigated by many authors (to name a few, see |31 ITTl IT2l fT51 HB] 
and references therein). One important step in this estimation procedure is the derivation of an 
asymptotic expansion of the voltage potential u in terms of the volume of the inhomogeneities. 
For an isolated inhomogeneity, the inhomogeneity dependent parameters in the first order of this 
expansion are called the Polarization Tensor (a constant 2x2 matrix) . The generalization of this 
tensor includes higher order terms dependent on the inhomogeneity via what is referred to as the 
Generalized Polarization Tensor (GPT) (see (3j [2j [5]) and was generalized to the case of linear 
elasticity (see [5]). 

These tensors appear in various contexts. In homogenization theory, they provide asymptotic 
expansions for dilute composites (see |171l5]). Very recently, GPT's have been used to reconstruct 
shape information from boundary measurements [7] and to construct approximate cloaking devices 
(Sj. In these last two papers, the authors consider the following conductivity transmission problem: 



1. Introduction 



div(7Vw) = in Q 

u = <p on d£l. 




<fiv((l D 7i + (l-l D )7b)Vu)=0 inR 2 , 



as \x\ — » 00, 
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where ji > and 70 > are constants, ljj is the characteristic function of a bounded domain 
D C M 2 with a Lipschitz boundary and a finite number of connected components, and H is a given 
harmonic function. In the absence of the inclusion D, u — H . Then, the far-field perturbation of 
the voltage potential created by D is given by 

t^d a T(x)M a0 d p H(x o ) as \x\ -> 00, 



(1.1) 



u(x) — H(x) 



00 

E 

l«l,|/3|=l 



where xq is the center of mass of D 7 and T is the fundamental solution of the Laplacian, 

1 



T(x) 



2tt 



In led 



The parameters a = (a\,a2) and /3 = (/3i,/?2) are multi-indices, and |a| = a\ + a.%. The real 
valued tensors M a p, which depend on D and on the conductivity contrast 71/70 are the so-called 
GPT. In bounded domains, expansions similar to (1.1 1 are also available, where T is replaced by a 



suitable Green's Function (see [3]). To fix ideas, let us assume xq is the origin. As it is observed in 
[3], the expansion ( 1.1 1 can be simplified when expressed in terms of harmonic polynomials which 
are real or imaginary parts of (x± + 1x2)™ where x = (xi,x^). Introducing 

a n (x) :— r" cos(n9) — 9? ((xi + 1x2)™) , b n (x) := r n sm(n6) — 3 {{x\ + ixzf 1 ) , 

if H can be written in the form 



00 



H(x) = H(0) + a n (H)a n (x) + /3 n (H)b n (x), 

n=l 

then it is shown in |3] that the expansion ( ] 1 . 1 [ ) can be written under the form 

u{x) — H(x) 



E 2 ^£L E ( M ( a "^) an(H) + M (a m b n ) f3 n (H)) 

m — l ' n—l 



Ei>m( x ) 
2irm\x\ 2m ^ 

m=l 1 1 n—l 



(M (b m a n ) a n (H) + M (b m b n ) (i n {H)) , as |x| 



00. 



The authors of [6J refer to M(P, Q), where P and Q are any of the harmonic polynomials a n or 
b n , as the contracted GPT. Such an expansion is particularly relevant if the inclusion is observed 
on a sphere Sr (of radius R) centered at the origin. In such a case, it is well known that a n and b n 
form a basis of L 2 (Sr) - this is the Fourier series decomposition. The tensor coefficients M (P, Q) 
are given by the formula 



(1.2) 
where 

(1.3) 



M(P,Q) = [ P(x)4> Q (x)da(x), 

JdD 



solves the transmission problem 



in D and M 2 \ D, 



on dD, 



dv 



71 _ 

7o dv 



dQ 



dv 
as I a; 



on dD, 



00, 



where v is the normal pointing outside of D, and where the subscript + (resp. — ) refers to the 
limit taken along the normal from the outside (resp. the inside) of D. 

The subject of this paper is to detail a simple algorithm to compute approximate values of the 
contracted GPT M(P, Q). We focus on the case when 71 is constant, and D is a simply connected 
inclusion with a C 1 boundary. In Section [2J we show that this problem can be written in a weak 
form as a system of boundary integral equations on dD 2 := dD x dD. Such an approach is well- 
known amongt numerical analysts interested in integral equations. We use an integration by parts 
formula that can be found in [141 118j . In Section [3] in the spirit of mixed finite-elements methods, 
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we decouple the boundary voltage potential from the normal flux and approximate them in the 
following way: 

4> Q \_ £ span (ai, ajvi.&tfi)) 

f d d d 9 , 

£ span — ai,— h,...,— a N2 ,— b N , 
\ ov Ov Ov ov 



dv 

for a given a fixed N\ and iV2, greater than or equal to the order of Q. We verify that the 
corresponding linear system is invertible. There are two advantages of this approach. The first 
is that since a n and b n are explicit polynomials, their gradient can also be computed explicitly. 
The precision with which the normal is calculated depends on the parameterization of D, which 
is independent of Ni and N 2 . The second is that the problem then involves the computation of 
4(A r 1 + N 2 ) 2 integrals on dD 2 , and then the inversion of a 4(iVi + N 2 ) 2 matrix. In other words, 
this decouples the order of the contracted GPT from the discretization of the boundary D. In 
Section [4] we discuss the results obtained for disks and ellipses, for which the exact values of the 
contracted GPT are known. Finally, in Section [5], we present a Matlab routine (named Myriapole 
on the Matlab Central exchange repository) that implements this approach. 

2. Weak boundary integral formulation of the transmission problem 

In this section, we derive the following variational formulation satisfied by the transmission 
problem. 



Proposition 2.1. Let <f>^ be the solution of the transmission problem (1.3). Introducing the 
notations 



Q ._ 



• r \dD - v-/> q v 

the pair (vP,v^) satisfies for any tp £ H 1 ^ 2 (dD) , 



7o 



(fc + 1) / — (y)j^(x)r(x-y)do-(x)do-(y)+2k / v »(y)^ x )—^ °L da(x)da(y) 

Jg D 2 OT OT Jg D 2 OV x 

= \l ~(x)<t>(x)d<r(x)- [ f-(y)^) dT{ ^ y) da(x)da(y), 
2 JdD ov J 9D 2 dv dv x 

and for all ip £ H^ 2 (dD), 

(k + 1) / vQ(y)ij(x)r(x - y) da(x)da(y) - 2 / v9(y)1>{x) 9T ^ - V) da(x)da(y) 
JdD 2 JdD 2 ° v y 

(2.1) = -/ ^-{y)^{x)T{x-y)da{x)da{y). 



dD< 



dv 



Proof. In order to solve the transmission problem of (1.3 1, we separate the interior potential from 
the exterior by defining the following: 

(2.2) u Q {x) := 4> Q {x) for x e D, and Q cxt (x) := cj) Q (x) for x e R 2 \ D. 
We then introduce the gradient fields as 

(2.3) v« int (x) := VuQ(x), and v« oxt (x) := V0 Q ext (a;), 

and set = v^int ■ v on dD. A representation of the solution using Green's identities leads to 
the expression of the interior solution and its gradient in terms of first and double layer potentials. 
For all x £ D, we have 

v9(x) = / v9{y) dT{x ~ V) - v<*(y)r(x - y) da(y), 

JdD av y 

v Q int (x) = [ u q (y)V e ( dT{ * ~ V) )- v<3(y)V x T(x - y) da(y), 

JdD \ av V J 
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As x — > dD from the interior, using the jump conditions satisfied by single and double layer 
potentials [TS], we obtain for x £ dD, 



(2.4) 
(2.5) 



u Q (x) 
v Q (x) 



;.)D 



^(y) dT{ * V) - vQ(y)T(x - y) da(y), 



u Q (y) .. .. 
dD dv x dv y 



dv y 



dv~ 



da{y). 



Similarly, for all x e M 2 \ D we have 



Q {y) dT(x^ y) _ {vQ ^ {y) _ v)v{x _ y)Myl 



dD 

dT(x - y) 



y ext(z/)V ;E 



i)D 



(v Q ox t(y) • v) V x T(x-y)da(y). 



And passing to the limit as x -> dB for the exterior, we obtain 



(2.6) 
(2.7) 



1 



Q cxtM = -/ <t> Q e*{v) ^ - - (v Q cxt(2/) ■ u)T(x-y)da(y), 

JdD av y 



t{x) ■ V 



du, 

Q d 2 T(x-y) 
dD av x dvy 



r Q M ■ v) dT t. V) My), 



dv x 



We subtract §2J§ from and use ([O) and ([2^2} to obtain for x e dD 
(2.8) 2 



u Q (y) dTi * V) da(y) - (1 + k) I v^(y)T(x - y) da(y) 



i)D 



dv„ 



dD 

dQ 



(2.9) 



2k 



OD 



dD dv 

v Q {y f^2-y> da(y)-(l + k) [ U Q(y) 



Similarly, from (2.5 1 and ( |2.7| we have, for x £ dD 

dT(x-y) 



dv x 



(y)T(x-y) da(y) 
d 2 T(x - y) 



OD 



dv x dv y 



da{y) 



ldQ 
2 dv 



(x) 



We test (^8j against (j) € H l / 2 (dD) and obtain 
2k 



da(x)da(y) + (1 + k) 



OD 2 



1 



dQ 



{x)(p(x) da(x) 



dQ , dT{x-y) 
dD dv dv x 



(y)jT 0) r - y) da{x)da{y) 
;, D * ot dr 

-K-{y)n x ) a da(x)da(y). 

dD 2 du dv x 



2 JdD dv 

We have used the following integration by parts formula proved in [141 118| 



dD 2 



u Q {y)(j)(x 



d 2 T(x - y) 
dv x dv y 



da(x)da(y) = 



du d(fi 
dD 2 dr dr 



Finally, we test ( |2.9[ ) against ip e H 1 ^ 2 (dD) to obtain 
2 



u°(y)V(*)^ir -da{x)da(y) - (1 + jfe) / v Q (y)xjj(x)T(x - y) da{x)da(y) 



dD 2 



dv,. 



dD 2 



dQ 



(y)tp(x)T(x - y) da(x)da(y). 



dD 2 



which concludes the proof. 



□ 



To conclude this section, we provide two equivalent formulations for the tensor M(P,Q) in 
terms of the unknowns and v®. 
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Proposition 2.2. The contracted GPT defined in ( |1.2[ ) is also given by 

M(P, Q) = (k - 1) J P{x) (j£(x) + (k- l)«°(o:)) da(x), 

and 

f BP 

M(P,Q) = {k-l) / —(x)(Q + (k-l) u Q{x))da(x). 



JdD 9iy 

Proof. The first formula is proved in [?J Lemma 3.3]. The latter is obtained after an integration 
by parts of the former. □ 

3. Approximate resolution method 



To solve the system given in (2.1|, we restrict it to a finite-dimensional space. More precisely 
given Ni and N 2 two positive integers, we set n p := 2N% and n q := 2N 2 and define 

q.ti when i < n p is odd, 



Pi := 
and 

:= 



b± when i < n p is even, 



Voi+i • v when i < n q is odd, 
V6i • when i < n q is even, 



and we consider (2.1 1 as a system posed on 

"H = span (>!,..., P„ p ) x (Q 1} . . . , Q n J 

We write 

2JVi 2N 2 

u Q = C « P «' and yQ = d nQn, 



and thus (2.1 1 becomes a linear system of the form AiX — B, where X is a (n p + n q ) x n q matrix, 
where for q < k < N 2 , the fc-th column contains the coefficients the of (u ak , v°" k ) and the N 2 + fc-th 
column the coefficients of (u bk , v ak ) . For 1 < i,j < n p , 

f dP dP 

Mij = (fc + ljT'ij, with = / — ^(y)— ^(x)r(x - y) da(x)da(y). 
For 1 < i < n p and n p + 1 < j < n p + n q 

dT(x-y) 



Mij = 2fcMj, with Mj = / Qj- n Jy)Pi(x) — h. — da(x)da(y). 

JdD" 1 av x 

For n p + 1 < i < n p + n q , and 1 < j < rz p , 

Af w = -2 / Q J _ B> (g)P j (y) ar ^ da(x)da(y) = -2A/J,,. 
For n p + 1 < i, j < n p + n qi 

Mij =(k+ l)Qi.j, with Qi.j = / Qj_ n Jx)Qi- n (y)T(x - y) da(x)da(y). 

JdD 2 

The right-hand-side term B is a (n p + n 9 ) x n q matrix, given when 1 < i < n p by 



Bid = o / Qj(a;)Pi(x)do-(x) - / Qj(y)Pi(x) da(x)da(y), 

z JdD JdD 2 av v 



and when n p + 1 < i < n p + n q . 



Bi,j=- / Qj(y)Qj-n p (x)T(x - y) da(x)da{y). 

JdD 2 

Note that both V and Q are symmetric, negative definite provided \x — y\ < k < 1, that is, 
provided the perimeter of £> is smaller than 2k. This assumption is not restrictive: we can scale 
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D to be sufficiently small for the computations. In that case, since we established that the matrix 
A4 is of the form 

" (k + l)V 2kAf 

M = 

-2N T (k + l)Q 

it has as a positive determinant and therefore this linear-system is well-posed. For any i and 
j such that max(z,j) < min(n p ,n 9 ), the contracted GPT coefficients are then computed using 
Proposition |2.2| namely 



Mil',. IV. = (A - 1) J Pi(x) [ Q 3 {x) + [k - 1) ^ X ro , nj)+j Q m (x) ) do-fcc), 



or equivalently by 

M(P i} P s ) = (*-!) Q<(x) f + (A - 1) ]T X mj P m (x) J da(x) 



m— 1 



For a fixed pair increasing the number n p and n q improves the accuracy of the method, 

to a point. Note that P, Q and T are explicitly given in terms of x and v{x). To compute the 
integrals which appear in Ai, B and finally for M, we mostly need to choose a parameterization of 
3D, and to compute its normal. If dD is given by an analytic formula, these integrals can all be 
computed using symbolic computation software. In this sense, the method is partially algebraic. 
For non-explicit curves dD, the only non-routine step is the evaluation of the double integrals 
involving T(x — y) and W(x — y) ■ f(x) on the line-segments (or on the arcs) where x — y cancels 
(or is smaller than a threshold). However, approximations for these integrals can be computed 
"by hand," thanks to the explicit form of the integrands. 

4. Numerical results for disks and ellipses 

In order to test our method, we consider geometries for which exact formulae for the contracted 
GPT are known, that is, disks and ellipses. If the inclusion is a disk, the contracted GPT is 
diagonal, and given by 

(4.1) M(P 2i . 1 ,P 2j . 1 ) = M(P 2i ,P 2j ) = Sij-j^ O^j . 

The formula for ellipses is slightly more involved, but also explicit (see [5] Prop. 4.7]). In this 
case, an algebraic approach (as described above) would provide exact results. Thus we use these 
shapes as a first benchmark for more complex geometries, for which the analytic formulae do not 
exist or are unknown. We compute the integrals appearing in Section [3] using a midpoint rule 
and fixing N\ = N 2 . In a first experiment, for a perfect disk of radius .5, we study the influence 
of the resolution of the boundary discretization and the number of polynomials used in the basis 
to compute the contracted GPT up to order 4. Figure |4.1| shows a plot of the result. The value 
represented is the maximal relative absolute difference between the approximate tensor and the 
exact tensor. For a tensor M' approximating M, the relative error is given by 



(4.2) e = max 



|K-- M ttl 

i,j max \Mki\ 

fc,Z>min(i,j) 

As i and/or j increases, max |Mm| decreases: this scaling factor compensates the variation 

fc,i>min(i,j) 

of tensor's coefficient size. Figure |4.1| shows that the error behaves as we could expect. If the 
number of harmonic polymonial used is less than twice the order (in these results, n = 4), the 
results are incorrect. Note that the harmonic polynomials come in pairs as real and imaginary 
parts of (xi +ix2) m , thus it makes sense to expect accurate results for tensors of order n only if the 
maximal degree (m in this case) of the representation of u is greater than n (m > n). Once this 
threshold is passed, the error decreases with the number of discretization points independently of 
the number of polynomials used. 
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Log 1Q Maximum Relative Error 
Order = 4, Contrast = 1/3, R = .5 




5 10 15 20 25 30 

Number of Polynomials Representing u(x) 



Figure 4.1. Relative error for GPT up to order 4. 



With the same geometry and the same order, we use 9 harmonic polynomials, and test the 
behavior of this method as the contrast increases, for various discretizations. The result is shown 
in Figure |4.2| Higher contrast k > 1 leads to greater errors, but the error remains bounded when 
the contrast tends to zero. With the same geometry, and a contrast of 1/3, we now fix the number 




Figure 4.2. Error vs. conductivity ratio. The error is bounded for k < 1 and 
increases with k > 1. For a point-count of 1024 and k s» 10, the relative error is 
less than 1%. 

of polynomials used as a function of the order (2n + 1 polynomials, where n is the tensor order) 
and increase the order. The result is shown in Figure |4.3| The image on the right represents the 
same data as that on the left but with a different scale on the vertical axis to show more detail. 
When too few points are used to discretize the boundary, that is, 16 or 32 points, the left plot 
shows that the ill-conditioning of the linear system causes extremely divergent results. For the 
other cases, the results are consistent. For tensors up to order 28, a parameterization point-count 
of 256 points gives relative errors of less than 1 percent for this disk case. 

Moving away from a perfect disk, we now investigate how the algorithm fares with ellipses 
of high eccentricity. Again, fourth order tensors were calculated and the same error measures 
were used. The ellipse with semiaxis lengths of a = .005 and b = .5 gives vastly different results 
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5 10 15 20 25 30 5 10 15 20 25 30 



Tensor Order Tensor Order 

Figure 4.3. Results for higher order tensors for a disk of radius one and con- 
ductivity ratio equal to 1/3. 

compared to the results from the disk. The first difference to note is that with 32 points, the 
approximation is not even close to the exact tensor. 

The first plot in Figure |4.4| shows the expected general trend that for larger polynomial-count 
and boundary point-count, error is smaller. However, there is an interesting phenomenon for 
the case of high point-count; as the number of polynomials representing u and v increases, error 
attains a minimal value but then after a certain point begins to increase. This is most likely due 
to floating point errors in Matlab that round the high-degree polynomials off and then add this 
error across each boundary segment at points on the more pointed edges of the ellipse. 

The second plot in Figure |4.4| gives another clear example of the inaccuracies of the ellipse for 
32 boundary points, the increased error for k > 0, and the existence of an "optimal" polynomial- 
count for parameterizations of high point-count. Also, the difference between k < 1 and k > 1 is 
clear in this case, as can be seen in the second plot. 

Lastly, we considered the amount of time taken to calculate tensors for the unit disk. On a 
64-bit, 2.3 GHz machine, Matlab computed fourth order tensors according to the results given in 
Table [j] 



Table 1 . Calculation Time (seconds) : Timing results for a tensor of order 1 and 
10 each using two different polynomial-counts: one equal to the minimal required 
amount, and one equal to twice the minimal required amount. 



Tensor Order 


Number of Polynomials 


16 Points 


256 Points 


1024 Points 


1 


3 


.0149 


.0340 


.8890 




6 


.0188 


.0467 


1.1463 


10 


21 


.0308 


.1063 


2.7794 




42 


.3055 


.4712 


5.4148 



5. A Matlab package with graphical user interface 

To extend the algorithm, we have developed a Matlab graphical user interface that calculates 
contracted GPT's of any shape. The interface provides an efficient and user-friendly implemen- 
tation of the algorithm described above and allows the user to easily change parameters used in 
the calculation. In this context, it also serves as a useful tool for further investigation into the 
accuracy of the computation for disks and ellipses. It is available from the Matlab Central File 
Exchange server under the name 'Myriapole'. 
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Number of Polynomials Representing u(x) 

FIGURE 4.4. Results for ellipse of considerable eccentricity with semiaxes a = .01 
and 6 = 1. The general region where the relative error is below 10 percent is where 
polynomial-count is greater than 10 and point-count is greater than 200. 




Figure 5.1. The layout of the graphical user interface Myriapole. 



The principal parameters of contracted GPT's are its order and conductivity contrast. Order 
in this situation refers to the highest power of x\ + ix 2 appearing in the set of polynomials P 
used for the calculation of M(Pi, Pj). Therefore, a contracted tensor of order n will be a matrix 
of size 2n x 2n, where each odd row and column correspond to P m = $t((xi + ix 2 ) m ) while each 
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even row and column correspond to P m = $s((x-l + ix 2 ) m ). The tensor order and contrast, along 
with the boundary shape and size, are referred to as tensor properties and are separate from the 
other "approximation" parameters in the graphical user interface (GUI), namely the number of 
polynomials used to represent u or the number of points in the parameterization of the boundary; 
the tensor properties change the values of the exact tensor while the approximation parameters 
change the accuracy of the approximating tensor. 

The panels of the GUI separate its different functions. For disks and ellipses, one can use the 
interface to calculate both exact and approximate tensors. Exact tensors are calculated using 
analytical formulas while the approximate ones use the algorithm described in this paper. This 
is useful for testing the effects the input parameters have on the accuracy, and in fact a separate 
panel in the interface can be used to easily assess the errors of an approximately calculated 
GPT. For any exact tensor M and approximating tensor M', the included measures viewable in 
the GUI are the absolute difference \M' — M\, relative absolute difference and LP matrix norms 



Beyond the capability to test the algorithm using disks and ellipses, the GUI also allows the 
user to import and calculate the GPT of simple shapes described in a bitmap image file. The 
easiest way to calculate tensors for shapes in this manner is to draw the filled shape in black on 
some white background in a graphic utility and then import this image into the interface. A quick 
algorithm parameterizes the shape and thus allows the user to calculate its approximate tensor. 



We have introduced a simple method to compute contracted Generalized Polarization tensors, 
together with a graphical user interface to test it. We verified on benchmark cases, disks and 
ellipses, that the resolution method was accurate, and was quite robust in the case of ellipses with 
very large aspect ratio. This code is available to anyone, and can be downloaded from the Matlab 
Central file exchange repository. Paired with other routines that complete the asymptotic expan- 
sions of domains marked with inhomogeneities of contrasting conductivity, this Graphical User 
Interface or the underlying code can function as a useful tool in implementations topological op- 
timization methods, inverse conductivity tomography algorithms, calculations involving effective 
conductivity of dilute composites, and others. It is possible that this approach can be extended 
to compute anisotropic generalized polarization tensors. We hope that this semi-algebraic ten- 
sor calculation method will inspire similar numerical methods for related transmission and layer 
potential problems. 
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